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This work presents new parallelizable numerical schemes for the integration of Dissipative Particle 
Dynamics with Energy conservation (DPDE). So far, no numerical scheme introduced in the literature 
is able to correctly preserve the energy over long times and give rise to small errors on average 
properties for moderately small timesteps, while being straightforwardly parallelizable. We present 
in this article two new methods, both straightforwardly parallelizable, allowing to correctly preserve 
the total energy of the system. We illustrate the accuracy and performance of these new schemes 
both on equilibrium and nonequilibrium parallel simulations. 


I. INTRODUCTION 

Molecular Dynamics (MD) nowadays easily al¬ 
lows the simulation of systems composed of several 
millions of atoms for molecular systems. This is 
however still insufficient when interesting phenom¬ 
ena occur on space and time scales spanning sev¬ 
eral orders of magnitude. One such example is the 
simulation of shock waves in molecular materials, 
where the resolution of atomic vibrations limits the 
timestep to fractions of picoseconds when covalent 
bonds are present. Coarse-grained models, which 
reduce the complexity of MD simulations with full 
atomistic details by representing some degrees of 
freedom in an average manner, are therefore of pri¬ 
mary interest in this context, they allow to extend 
the spatial scale by reducing the number of degrees 
of freedom, and also the temporal scales by increas¬ 
ing the admissible timesteps since mesoparticles in¬ 
teract via smoother potentials. 

The Dissipative Particle Dynamics (DPD)i is a 
particle-based coarse-grained model in which atoms, 
molecules or even groups of molecules are repre¬ 
sented by a single mesoscale DPD particle. The time 
evolution of the mesoscale particles is governed by a 
stochastic differential equation. Dissipative and ran¬ 
dom forces allow to take into account some effect of 
the missing degrees of freedom. DPD was put on a 
firm theoretical ground in Ref-d, and was applied to 
study the properties of various systems^ - —. 

However, DPD intrinsically is an equilibrium 
model, with a prescribed temperature. The temper¬ 
ature is fixed through a fluctuation/dissipation re¬ 
lation between the magnitudes of the friction and 
the fluctuation parameters. DPD cannot be used as 
such to study nonequilibrium systems. It should 
be replaced by a dynamics where the fluctuation- 
dissipation relation is based on variables which 
evolve in time. DPD with conserved energy (DPDE) 
is such a model^ 7 . 

In the DPDE framework, mesoparticles have an 


additional degree of freedom, namely an internal 
energy, which accounts for the energy of the miss¬ 
ing degrees of freedom. The dynamics on the inter¬ 
nal energies is constructed in order for the total en¬ 
ergy of the system to remain constant. DPDE was 
for instance used to simulate shock and detonation 
waves& - — . 

However, the efficient numerical integration of 
DPDE still requires some effort. Numerous effi¬ 
cient schemes were developed for DPD 11 . - — (see in 
particular Ref. [Til for a careful comparison). How¬ 
ever, their adaption to DPDE leads to numerical 
schemes for which errors on average properties may 
be large even for timesteps standardly used to inte¬ 
grate Hamiltonian dynamics^. This may therefore 
require the use of extremely small timesteps 17,18 . 

Among the currently known schemes for DPDE, 
one of them turns out to enjoy very nice properties 
in terms of energy conservation. It is a splitting 
scheme inspired from Shardlow's splitting scheme 
for DPD*™, and adapted to the DPDE framework^. 
The resulting scheme, called SSA in this work 
(Shardlow Splitting Algorithm), is so far considered 
as the reference scheme for the numerical integra¬ 
tion of DPDE^l. The main issue with SSA is that par¬ 
ticle pairs have to be updated sequentially. The par¬ 
allelization of SSA is therefore not an easy tasl<—■. In 
addition, SSA is not threadable, and thus incompati¬ 
ble with the architectures of future supercomputers. 

The aim of our work is to propose and test new 
schemes to integrate DPDE. We pay a particular at¬ 
tention to (i) the preservation of energy, (ii) the size 
of the errors on average properties arising from the 
use of finite timesteps, and (iii) the possibility to eas¬ 
ily parallelize the method. More precisely, our aim is 
to construct easily parallelizable schemes for which 
the systematic errors (i.e. the remaining errors in the 
limit of infinite sampling precision) are as small as 
possible. 

Let us also emphasize that the schemes we de¬ 
velop for DPDE may be of direct interest for 
other dynamics with similar structures such as the 
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Smoothed Dissipative Particle Dynamics (SDPDj^i. 
This stochastic dynamics also preserves the energy 
of the system. It can also be seen as the superpo¬ 
sition of two elementary energy preserving dynam¬ 
ics, a conservative part and a fluctuation/dissipation 
part. It is therefore not a surprise that the numer¬ 
ical methods we develop here can be readily ap¬ 
plied for large scale simulations of SDPD on parallel 
architectures 22 . 

This article is organized as follows. We first briefly 
recall the DPDE model and its properties in Sec¬ 
tion El We next present the numerical schemes we 
consider in Section um some reference schemes in 
the literature are recalled in Section IIII A1 while 
new schemes are proposed in Section fill Bl We then 
perform sequential and parallel equilibrium simula¬ 
tions to evaluate the quality of the energy conserva¬ 
tion and the errors on average properties as a func¬ 
tion of the timestep (see Section HVll . Finally, in order 
to test the dynamical relevance of the newly intro¬ 
duced schemes, we present parallel simulations of 
shock waves in Section IVl 


II. PRESENTATION OF THE MODEL 


set of equations: 


dq ilt = v i>t dt, 

N 

dpi,t = -V qi U{q t )dt- Y yij,tX( r ij,t)Vij, t dt 

j=bj*i 

N 

+ Y crJx(rij,t)dW ij>t , 

M.M ( 1 ) 

N i 2 \ 

N 

- Y a ^x( r ij,t) v ij,t ■ dw ij,t > 

where (^ij,t)\<i<j<N are standard d-dimensional 
Brownian motions, and W t j t = -Wj l t for i > j. In the 
above equations, we have introduced the relative ve¬ 
locity Vijj = v lt - Vj F f between particles i and j and 
their distance r lht = \\q it - qj t \\. The function x( r ) is 
a cut-off function, limiting the range of the random 
and dissipative interactions. We consider here 



X(r) 



if r < r cut , 


otherwise. 


( 2 ) 


A. The DPDE equations 


We consider a system of N identical particles in di¬ 
mension d, with positions q\ and momenta p;. Their 
velocities are Vj = fk, where m is the mass of the par¬ 
ticles (masses are taken to be identical to simplify 
the presentation; the extension to particles with dif¬ 
ferent masses is straightforward). We rely on the 
DPDE models 7 although alternative deterministic 
models have also been proposed^. In the DPDE 
framework, each particle represents a (group of) 
molecule(s), p, and p, being the position and momen¬ 
tum of the center of mass of this collection of atoms. 
These external degrees of freedom are coupled with 
some internal energies £; describing in an average 
manner the missing degrees of freedom lost in the 
coarse-graining process. Each pair of particles in¬ 
teracts through (i) conservative forces deriving from 
a potential energy function U(q), (ii) friction forces 
proportional to the relative velocity between the two 
particles and (iii) random fluctuation forces. In this 
article, we consider a variant of the original DPDE 
model introduced in Ref. @L where friction and ran¬ 
dom forces are not projected along the line joining 
two particles. However, the numerical schemes pre¬ 
sented in Section Hill can be straightforwardly modi¬ 
fied for the original DPDE model. 

The time evolution of the configuration ( qi,pi,£i ) 
of the ith DPDE particle is given by the following 


Note that, by construction, the total momentum 
Yi =1 Pi and the total energy 


N 2 N 

£(q,p,e) = U(q) + Y— + Y e ' 

l —' m l —■ 

i =l i=i 


are preserved. 

The parameters and a respectively control the 
friction and fluctuation strengths. They are chosen 
in order to ensure the correct statistical behavior of 
the system. An important point is that the invariant 
measure should be consistent with the preservation 
of total energy and momentum. More precisely, 
and a are chosen in order for measures 


. N \ 

f[£{q,p,e)}g Yp i exp^<S(£)jdgdpd 

i=i i 


( 3 ) 


to be invariant under the dynamics CQ for arbitrary 
functions / and g. In this expression, 


N 

■S(z) = Y s ^ 
1=1 


is the total internal entropy of the system, s(e;) being 
the internal entropy of the ith particle. The internal 
entropy function satisfies 
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where T(e) is the internal temperature implicitly de¬ 
fined via the relation 


e = 


rT(e) 

Jo ' 


C v (6)d6. 


( 4 ) 


Here, C v (0 ) is the internal heat capacity at constant 
volume. Therefore, the internal heat capacity fully 
determines the internal energies and the internal en¬ 
tropies. The relation is the internal equation of 
state of the DPDE particle. 

In order for measures of the form ([3j to be invari¬ 
ant by the dynamics CD, the parameters y t j it and a 
should satisfy the following fluctuation/dissipation 
relatior 6,7 - 


m 


o 


1 1 

4 \7) + 7 ) 


( 5 ) 


where T; = T(£,) is the internal temperature of 
the zth particle. In order to highlight the similar¬ 
ity with the standard fluctuation/dissipation rela¬ 
tion for Tangevin dynamics (see for instance Sec¬ 
tion 2.2.3 in Ref. ~24l) . we find it convenient to rewrite 
the definition of y,j and a as 


where Zc„ d„ is a normalization constant. Under this 
assumption, average properties can be computed as 


< A > = J AA}lEo ’ P ° = fJi’SoT J Q A (^Pr^x)du 


( 7 ) 


where ( q T ,p r ,e T ) is the solution at time T of CD start¬ 
ing from a given configuration (qo,po,£o)- This ini¬ 
tial configuration is such that £(q 0 ,p 0 ,£ 0 ) = E 0 and 

Hill Pi,o = Pq- 

The measure Pe 0 ,p 0 can be interpreted as some mi- 
crocanonical measure. In order to obtain estimates 
of thermodynamic properties such as temperature, 
it is convenient to introduce a canonical equivalent 
of the measure Pe 0 ,p 0 • In fact, the measure Pe 0 ,o (cor- 
responding to a system with total momentum set 
to 0) should be equivalent in the limit of large sys¬ 
tems to the canonical probability measure 

dpp (q,p,e) = d q dp de, (8) 

where Zp is a normalization constant. The parame¬ 
ter is chosen in order for the average energy un¬ 
der (8| to be equal to E 0 : 


Pd 

7ii = 7 ^’ 



( 6 ) 


J £(q,p,£)dpp(q,p,e) = E 0 . 


with 

Pij = + Tjr 

/c B being Boltzmann's constant. This amounts to in¬ 
troducing a reference friction parameter y and a ref¬ 
erence temperature T re f, with associated inverse en¬ 
ergy &ef = I?h- 

B. Thermodynamic properties 

In order to compute the average value of a physi¬ 
cal observable A, we assume that the dynamics CD is 
ergodic for some invariant measure p of the form ([3j. 
Note indeed that there are no mathematical results 
ensuring the ergodicity of the stochastic dynam¬ 
ics CD since the fluctuation part may be degenerate 
(contrarily to, say, Langevin dynamics). Ergodicity is 
even not known to hold for DPD, except in the very 
specific case of one-dimensional systems^. 

Since the total energy and the total momentum 
are preserved by the dynamics, the only reasonable 
assumption is that ergodicity holds with respect to 
the probability measure 

PE 0 ,p 0 (dqdpde) 

= Z £ 0 1 ,P 0 ,5 (£:( ? ,p, £ )-£ 0 }<5j L N ip ._ Po |e' S(t ) 7 


This equivalence is very similar to the standard 
equivalence between microcanonical and canonical 
measures for systems described only by their posi¬ 
tions and momenta. 

We focus in the numerical tests we present on var¬ 
ious estimators of the temperature, each one involv¬ 
ing only one type of the three categories of degrees 
of freedom of the system: the positions p,-, the mo¬ 
menta p,- and the internal energies £;. The interest 
of these estimators is that they should all predict 
the same temperature (provided the system is suf¬ 
ficiently large, which, in our experience, is already 
the case for moderately large systems of several hun¬ 
dreds of particles). This therefore allows to assess 
the quality of the sampling for each group of degrees 
of freedom. 

The first estimator of the temperature is the stan¬ 
dard kinetic temperature 



where (-)p refers to averages with respect to the mea¬ 
sure pp introduced in ©, and N e $ represents the 
effective number of external degrees of freedom of 
the system. It is a priori equal to dN but since we 
fix the total momentum to Pq = 0, we reduce it to 
N e ff = d(N - 1). This correction is anyway unimpor¬ 
tant for sufficiently large systems. 
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The second estimator— depends only on the posi¬ 
tions of the particles: 


, (ll^lP), 

p0t h (AU)p 


( 10 ) 


Finally, the last estimator is determined by the inter¬ 
nal energies: 



A simple computation shows that Tpot - lint — Jkin - 
k B //I (the second equality is proved in Ref. \27\). Let 
us emphasize that the internal temperature is esti¬ 
mated by a harmonic average rather than an arith¬ 
metic one. 


C. Practical computation of average properties 

In order to estimate in practice average proper¬ 
ties by the ergodic limit {7J, the dynamics needs to 
be numerically integrated. We therefore introduce 
a timestep At > 0, and denote by ( q n ,p n ,c n ) an ap¬ 
proximation of the solution ( ^nAt’PnAt’^nAt) °f CD 
obtained by iterating a numerical scheme. Average 
properties are then estimated by simulations over 
N[ tel timesteps as 

< tq ter 

<A)-A A , Nte = — YA(q’',p n ,e"). (12) 

JViter n=1 

Standard results from the numerical analysis of 
stochastic differential equations (see for instance 
Section 2.3.1 in Ref. f24l) allow to quantify the errors 
produced by the approximation (1121) as follows: 

+ , 13 ) 

This equality highlights two sources of error in the 
computation of average properties: 

• a statistical error R AAt /yJN[ ter At, which arises 
from the finiteness of the number of iterations. 
Typically, a central limit theorem holds, so that 
R A ,At follows some Gaussian distribution with 
a variance (? AAt close to the variance of the 
time averages estimated with the underlying 
continuous dynamics. The important point is 
that the statistical error decreases as the in¬ 
verse of the square-root of the physical simu¬ 
lation time. 

• a systematic bias K A At’t, which is the residual 
error persisting in the limit of infinite sam¬ 
pling accuracy. This error arises from the use 
of finite timesteps. 


The statistical error only mildly depends on the nu¬ 
merical scheme at hand since the variance °l,At is 
at first order in At the variance of the continuous 
process— As in Ref. [ 20 L our focus in this work is 
therefore rather on the systematic bias, which may 
be quite different for various numerical schemes. 
Our aim is to construct schemes for which the sys¬ 
tematic bias is as small as possible. For a given maxi¬ 
mal number of step N iter , this allows to integrate the 
dynamics with larger timesteps, and hence reduce 
faster the statistical error. 


III. INTEGRATING DPDE 


The numerical schemes we consider in this section 
are based on splitting techniques as presented in the 
DPD context— ?,2S . The idea of splitting algorithms 
is to decompose the dynamics into several elemen¬ 
tary subdynamics, which are sequentially integrated. 
Splitting schemes are interesting when the elemen¬ 
tary dynamics are simple to numerically integrate. 

Although splitting schemes for stochastic dynam¬ 
ics are nowadays quite popular, some standard inte¬ 
gration schemes for DPD do not fall into this class. 
For instance, some integrators were obtained by an 
ad-hoc modification of the Verlet schema^ for the 
integration of Hamiltonian dynamics— For com¬ 
pleteness, we also consider in our numerical exper¬ 
iments the adaption of the Velocity-Verlet scheme to 
the DPDE context. This scheme is abbreviated SVV 
in the sequel for Stochastic Velocity-Verlet, and is 
made precise in Appendix lAl 

Let us also immediately emphasize that a compli¬ 
cation of DPDE compared to other stochastic dynam¬ 
ics such as DPD arises from the presence of two 
invariants of the dynamics, namely the total mo¬ 
mentum and the total energy. Ideally, numerical 
schemes should preserve these invariants over very 
long times, at least approximately (as the energy is 
approximately preserved by appropriate discretiza¬ 
tions of Hamiltonian dynamics). 

We first briefly present in Section fLII Al two split¬ 
ting schemes already known by the community: a 
simple splitting scheme based on a Euler-Maruyama 
discretization of the dissipative part of jT} (termed 
SEM in the sequel) and the adaption to the DPDE 
context of the splitting proposed by Shardlow in the 
DPD contexfcM? (termed SSA in the sequel). The 
simple splitting scheme is very easy to parallelize 
but leads to very large errors in the estimated proper¬ 
ties; while SSA is quite accurate but somewhat cum¬ 
bersome to parallelize— This motivates the intro¬ 
duction of two new straightforwardly parallelizable 
schemes in Sections lIII B ll and llll B 21 
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A. State of the Art 


The DPDE dynamics {TJ can be decomposed into 
a Hamiltonian evolution 


dq, t = dt, 
m 

K d Pi,t = -V qi U(q iit )dt, 
and a fluctuation/dissipation part 

N 


d Pi,t = - Y 7ij,tX(rij,t)Vijdt, 
i=hj*i 
N 

Y a ^X(r l] ,t)dW lht , 


de U = \ 


j=l,j*i 
N 


j=7j*i 

N 


Y yd-^ij.t ~ d—\x{ r ij,t)dt, 


(14) 


(15) 


- Y a ^X( r ij,t) v ij,f dW ij,t 


j=hj*i 


All splitting schemes first integrate the Hamiltonian 
part with a standard Velocity-Verlet discretization- 29 . 
The difference therefore solely arises from the subse¬ 
quent treatment of the fluctuation/dissipation part. 


1. Splitting Explicit-Euler: SEM 


The SEM scheme integrates the fluctua¬ 
tion/dissipation with a simple Euler-Maruyama 
discretization: 


N 


pr 1= Pi- Y 

j=l,j*i 

N _ 


+ Y 

j=\,j*i 

1 N i 2 \ 

cr‘^r+2 L mK i 2 -^b<d> A ' 

i=Ei*i ' ’ 


(16) 


1=1,] 

N 


\ Y 


i=hj*i 


where (here and in the sequel) (G,” )neiN,i<i<;'<N 
are identically and independently distributed stan¬ 
dard d-dimensional Gaussian random variables, and 
G"j = -G”.. Note that no particular care is taken to 
make sure that the energy is preserved by the dis¬ 
cretization of the fluctuation/dissipation part. Un¬ 
surprisingly, it turns out that the total energy drifts 
in time (see Section llV Bl) . 


As SVV, SEM turns out not to be a very accurate 
scheme. It can be shown to be of weak order one. 
Under appropriate ergodicity conditions, the bias on 
average properties is therefore of order At (;/ = 1 in 
HU). We refer to Appendix[B]for more precisions. 

The reason why we consider SVV and SEM in our 
numerical tests is that they can be efficiently paral¬ 
lelized in a straightforward manner, and are there¬ 
fore appealing alternatives to SSA for large scale sim¬ 
ulations on massively parallel architectures. 


2. Shardlow’s Splitting Algorithm: SSA 


SSA, in opposition to SEM, further decomposes 
the fluctuation/dissipation dynamics (1 1 51) into ele¬ 
mentary pairwise fluctuation/dissipation dynamics 
involving only two DPDE particles i and j. Tak¬ 
ing advantage of the energy and momentum con¬ 
servations, these elementary fluctuation/dissipation 
equations can in fact be rewritten as follows: 


d Pi,t = -Yij.tVij,tX{rij,t)dt + afyjx(rij,t)dWij iU 
d Pi,t = ~ d Pi,t’ 


„2 


d£j, f =~2 d 

d £j,t — dcif. 


2 \ 


Ek + p Jl 

2m 2m 


(17) 


As noted in Ref. [28| for DPD, the dynamics on the 
momenta can be analytically integrated when y t j is 
fixed (i.e for fixed ej), instead of using a modified 
BBK algorithm^ as done in Ref. ITil . For the exten¬ 
sion of SSA to DPDE we consider in this work, we 
rely on the strategy introduced in Ref. 0, where mo¬ 
menta are updated first (here with an analytical in¬ 
tegration at fixed y,y), and internal energies are then 
updated to ensure the overall energy conservation. 
Introducing 




2 7ij 

m 


X(ru)At, 
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the scheme we use to integrate the elementary pair¬ 
wise fluctuation/dissipation (1 1 71) reads 


pr +1 =p;S 


+ CM 


pr=Pr - 2 


+ a -1 



-M 

ml 

- e ~ 2 a "j,At j 

J 




ml 

-e“ 2a P- Af j 



r n 

G ij 


( 18 ) 


sequel (see SectionAs parallel simu¬ 
lations are performed using a spatial reparti¬ 
tion of the simulation box between the proces¬ 
sors, the bottleneck for the parallelization of 
SSA arises from particles located in different 
domains^. Therefore, the idea of the Hybrid 
scheme is to integrate the elementary fluctua¬ 
tion/dissipations interactions involving parti¬ 
cles located on the same processor by a pass of 
the SSA algorithm, while the remaining inter¬ 
actions are taken care of by a SER discretiza¬ 
tion. 


1. Splitting with Energy Reinjection: SER 


„n+l 


= *?" 2 AK ij' 


c n +l - C n A V n 

e i ~ £ i - 2 H’ 


The SER integration of the fluctua¬ 
tion/dissipation (1151) is performed in two steps. 
First, momenta are integrated using a simple 
Euler-Maruyama discretization as 


where we introduced the kinetic energy variation 

AK ij = ^ [(pr 1 ) 2 + (Py +1 > 2 - (PD 2 - (P") 2 ] 

in order to update the internal energies. Particle cou¬ 
ples which are sufficiently close to interact are then 
sequentially updated with the scheme d 1 81) . It is pre¬ 
cisely this sequentiality which prevents a simple par¬ 
allelization of the scheme. As a side remark, note 
that when linear scaling techniques are used (such 
as decomposing the system into cells of size r cut ), the 
order in which the particles are updated may change 
from one iteration to the other. 


B. New Numerical schemes for integrating 
DPDE 

The two new parallelizable schemes we propose 
rely on the splitting between Hamiltonian and fluc¬ 
tuation/dissipation parts mentioned at the begin¬ 
ning of Section fill A1 but with new strategies to dis¬ 
cretize the fluctuation/dissipation part: 

• The first scheme we present in Section fill B II 
is called “Splitting with Energy Reinjection" 
(termed SER in the sequel). Its discretization 
of the dissipative part {5} is similar to SEM 
but uses a global symmetric reinjection of the 
kinetic energy variation into the internal en¬ 
ergies, instead of directly discretizing the dy¬ 
namics of the internal energies. This allows to 
automatically preserve the total energy during 
the fluctuation/dissipation part. 


s P ?= b P ?, (i9) 


with 5p» = -yfjX^v^At + aJ^)Gf j y/Kt. The in¬ 
ternal energies e, are then updated in order to com¬ 
pensate for the energy variation during the update 
of the momenta. In order to implement this idea, 
we need to identify in the global kinetic energy vari¬ 
ation of each particle the contribution of every sin¬ 
gle pairwise interaction, in order to redistribute the 
associated elementary energy variations. In fact, a 
simple computation shows that 


AK" = 


(p| ,+1 ) 2 

2m 


(Pi 1 ) 2 

2m 


N 


L a a- 


where 



represents the contribution of particle j to the ki¬ 
netic energy variation of particle i. The internal ener¬ 
gies are then updated by reinjecting the elementary 
variations A jK" in a symmetric manner: 

, , (20) 


The second scheme is a mix between SSA and 
SER, and is therefore termed “Hybrid” in the 
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The discretization of equation <11711 can therefore be 
summarized as 


' N 

p'r l =p'h Y, rW r V<i At 

j=hj*i 

N 

+ L IxtyGij'fot- 

j=l,j*i 

N 1 

'?*■=<- E ^rHh^K- 4 '’"))- 

, j=bj*t 

( 21 ) 

It can be shown that this numerical scheme is con¬ 
sistent. In fact, it has a weak order 1, so that the 
bias on average properties is of order At. This can be 
checked by some straightforward but lengthy calcu¬ 
lus, see Appendix lB 21 

Let us emphasize that, by construction, SER au¬ 
tomatically ensures the exact conservation of the to¬ 
tal energy during the numerical integration of the 
fluctuation/dissipation part (II 5> . This very nice fea¬ 
ture was only enjoyed by SSA among the previously 
known schemes to integrate DPDE. On the other 
hand, in opposition to SSA, SER does no rely on a 
further splitting of the fluctuation/dissipation into 
elementary pairwise interactions, and is therefore 
straightforward to parallelize. However, the pres¬ 
ence of both dp”. and dp" in <12011 requires two sweeps 
on the system, a first one to compute dp"j and sum it 
into dp", and a second one to compute the products 
dp"j ■ dp". This seemingly requires to store all incre¬ 
ments dp"j (or at least the Gaussian increments G” ), 
which is somewhat prohibitive. This can however 
be avoided if the Gaussian increments can be ex¬ 
actly regenerated to the value they had on the first 
pass of the algorithm when the particle pairs are re¬ 
visited. The SARU pseudo-random generator^- pre¬ 
cisely allows such an easy recomputation since the 
Gaussian increments G”- are completely determined 
by the integers i,j and n. 


The variable G = (G^) (< ;<^ refers to the Gaussian 
increments used. for the discretization. We also 
denote by (q,p, e, G) the result of a SSA dis¬ 

cretization (Tb) of the elementary dynamics (II7> 
with timestep At. Note that only the components 

Pi,Pj,Ei,Ej of ( q,p,E) are changed by ® Af JI (q,p,E,G). 
Let us denote by A the set of particle couples where 
both particles are in the same processor, and A c 
the remaining couples. The elements of A are de¬ 
noted by (ii,ji), (i 2 >ji)> ■■■> where / is the 

number of elements in A. Finally, we denote by 
® Af ’ (q,p,E,G) the one-step iteration correspond¬ 

ing to a SER discretization of all the interactions in¬ 
volving the particle couples in A c . The expression 

of ® Af ybnd (q,p,e, G) is obtained from (q,p, e,G) 

and (q,p,£,G) as 


0 H y brid /~ ~ . /m - mSF.RW „ , 

At 


(q,p,E,G) = <i>lf < ’ o ®~(q,p,£,G), 


with 

®£(q,p,e,G) = $ 


SSA ,iiji 
At 


0...00 


SSA,i!,;i 

At 


(q,p,E,G), 


where o denotes the mathematical composition. 

The Hybrid scheme, in opposition to all the other 
schemes, does not depend only on the timestep but 
also depends on the space repartition of the simu¬ 
lation box between processors used for the simula¬ 
tions. This dependence is studied in Section llV El 


IV. EQUILIBRIUM SIMULATIONS 

In this section, we perform sequential and parallel 
equilibrium simulations using all the schemes pre¬ 
sented in Section IIIII However, the Hybrid scheme 
can only be tested in parallel simulations, otherwise 
it reduces to SSA. Therefore, no Hybrid scheme is 
used in the sequential simulations of Section IIV Bl 
and llVD II 

A. Description of the simulation 


2. A scheme between SSA and SER: Hybrid 

The Hybrid scheme can be seen as a blending 
of SER and SSA, where the elementary fluctua¬ 
tion/dissipation interactions involving particles on 
the same processor are integrated by a pass of the 
SSA algorithm, while the remaining interactions are 
integrated with the SER scheme. 

Let us describe more precisely this algorithm. We 

denote by 0 Af ybnd (^,p,£, G) the result of a Hybrid dis¬ 
cretization of CD. with timestep At, applied to the 
configuration ( q,p,E ) of the system, i.e 

(q n+1 ,p n+1 ,£ n+1 ) = ®*J hlid (q n ,p n ,£ n ,G n ). 


We consider a system of N particles, with pair¬ 
wise conservative interactions described by a shifted, 
splined and truncated Lennard-Jones potential: 




e sh fspi r r cut) 

if r < r cut 

0 

if r > r cut 


( 22 ) 


The parameters £ s h and / sp are chosen so that 
v(r cut ) = 0 and v'(r cut ) = 0. The total potential en¬ 
ergy of the system reads 

U(q)= Y v ( r ij)- < 23 ) 

0<i<j<N 
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In all the simulations presented below, the heat ca¬ 
pacity is supposed constant. Equation Q then sim¬ 
ply reduces to £,• = C v Tj. 

All results are given in reduced units and, un¬ 
less otherwise mentioned, the fluctuation strength 
is set to a = 4. The physical dimension d is set to 3 
with the particle number set to respectively N = 343, 
N = 2000 and N = 144,000 in Sections ITTBl IIV D 1 1 
and II V D~2l We fix C v = 50 for the sequential simula¬ 
tions of Section flV Bl and II V D H and C v = 10 for the 
parallel simulations of Section llV D 21 

Initial conditions are obtained in Sections II V Bl 
and IIV D II by melting a simple cubic crystal using 
a SSA scheme for t = 100 with a timestep At = 0.001. 
The initial condition for the parallel simulations of 
Section IIV D 21 is obtained by melting a cubic face 
centered crystal using a SER scheme for t = 45 with 
a timestep At = 0.001. The density is p = 0.8095 in 
all cases except in Section llV Bl where p - 0.5787. 


Figure |T] shows the behavior of the time- 
dependent energy drift averaged over trajectories, as 
a function of time, for various schemes and At = 
0.006. Note that no scheme manages to keep the 



B. Energy drifts 


The first task of an appropriate numerical scheme 
for DPDE is to preserve (at least approximately) the 
invariants of the dynamics, namely the energy and 
the total momentum. While the total momentum 
conservation is easily ensured, the energy conser¬ 
vation on the other hand is much more demand¬ 
ing. In fact, the numerical schemes we consider 
are obtained by a composition of a Verlet scheme, 
which approximately preserves the energy, and a dis¬ 
cretization of the fluctuation/dissipation dynamics, 
which may or may not preserve the energy. The in¬ 
teraction between these two integrators may lead to 
drifts in the total energy, even when the integration 
of the fluctuation/dissipation part exactly preserves 
the energy, as already observed in Ref. 1201. 

In order to precisely quantify the possible energy 
drifts and determine in particular the influence of 
the timestep At on the drift rate, we compute the av¬ 
erage evolution of the energy in time by performing 
N( ra j = 10 4 trajectories of time fy = 10 each. The ini¬ 
tial conditions of each trajectory are sampled accord¬ 
ing to the canonical measure ftp, obtained by sam¬ 
pling independently the internal energies accord¬ 
ing to the measure Z“ 1 e _ ^ £+S ^^, and the positions 
and momenta according to the canonical measure 
Z-y^'Pl We denote by (A t ) the value of an ob¬ 
servable A at time t, obtained by averaging over all 
the possible trajectories. In practice, (A nAt ) is ap¬ 
proximated by 


N. 


traj 


A n 

A£N traj 


N, ra ; / ' 

tra ) m =1 


A(x n 


(24) 


where x m,n is the n-th configuration of the m-th tra¬ 
jectory. 


Figure 1. Time-dependent energy drift averaged over tra¬ 
jectories, as a function of time when At = 0.006. 

total energy constant. Linear drifts in time are ob¬ 
served in all cases, SVV being substantially worse 
than all the other schemes in terms of energy con¬ 
servation. We study more precisely the drift rates in 
the remainder of this section, in order to determine 
the maximal timesteps which can be used in practice 
before average properties are corrupted. 


1. Quantification of the total energy drifts 

In order to decide whether the drift is acceptable, 
we quantify in this section the rate of variation of 
the energy as a function of the timestep At. Figure [I] 
suggests a linear fit of the time-dependent average 
energies as 


(£">Af = E o( l +a At nAt), (25) 

which allows to identify the relative energy drift rate 
a At- Figure [2] presents the dependence of a& t as a 
function of the timestep. 

The SER scheme has energy drift rates of the same 
order than SSA, whereas the drift rates for SVV 
are substantially larger in absolute value. In addi¬ 
tion, we notice that the drift rates of SSA and SER 
schemes remain lower than 5 x 10 -6 , except for SER 
at At > 0.007. This means that the energy vari¬ 
ation would be less than 5% for simulation times 
lower than f max = 10 4 . Such simulation times may be 
long enough to estimate equilibrium properties, but 
they may prove somewhat short to estimate trans¬ 
port properties. 

Moreover, the results of Figure [2] suggest that the 
drift rate a has a polynomial behavior with respect 
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Figure 2. Average energy drift rate a&t> as a function of the 
timestep. A fit a^f = KAt K is superimposed to the data. 


to At : 



At 


Figure 3. Drift rates per unit time of the kinetic (top), po¬ 
tential (middle) and internal (bottom) energies. 


ctA t — KAt K . 

A least-square fit in a log-log diagram gives k — 2 
for SVV, and k — 4 for the other schemes. This fast 
increase of the drift rate places a severe limitation 
on the use of larger timesteps in DPDE simulations. 

2. Drifts of individual energy components 

We studied in Section IIV B II total energy drifts. 
Let us now have a more precise look at the drifts 
of individual energy components (kinetic, potential 
and internal), for which the drift rates as a function 
of the timestep are plotted in Figure [3] As expected, 
the drift rates for SVV are large for all energy com¬ 
ponents, and small for all components for SSA. SER 
on the other hand exhibits a non-trivial behavior: its 
drift rates are quite large for individual components 
of the energy, although these large drifts compensate 
each other in order for the total energy to drift only 
slowly. In fact, we even observed in some situations 
that the drift of some energy components was not 
linear as a function of time. The very notion of drift 
rate for SER is therefore dubious (see Ref. [32j for fur¬ 
ther precisions on this specific issue). 

C. Projected dynamics to ensure stability 

In order to remove the energy drift, a natural ap¬ 
proach, already suggested in Ref. [ 20 !, is to enforce 
the energy conservation by an appropriate projec¬ 
tion on a surface of constant energy. There are sev¬ 
eral possible projection procedures. 

One possible option is to resort to some Lagrange 
multiplier. A configuration x = (q,p,e) such that 


£(x) E 0 is replaced by x - AV£(x), with A chosen 
such that 

£^q - A V q £(x),p - AV p £(x),e - AV t £(x) j = E 0 , 
or more explicitly 

£{q-XVU(q),p-—p,E-Al) = E 0 , (26) 

where 1 is the N-dimensional vector whose compo¬ 
nents are all equal to 1. However, the numerical 
computation of the parameter A satisfying <12611 is 
a computationally expensive task typically requir¬ 
ing several iterations of a Newton procedure, and 
hence several evaluations of the energy per timestep. 
Note also that the total momentum is not preserved, 
so that this extra conservation law should be subse¬ 
quently enforced in an appropriate manner. 

It is therefore much more convenient from a prac¬ 
tical viewpoint to play on the internal energies only 
to adjust the total energy - as was also already sug¬ 
gested in Ref. [2ol . More precisely, at the end of one 
iteration of the numerical scheme, the resulting in¬ 
ternal energies are replaced by e” +1 - (£(x”) - E 0 )/N, 
where x n is the unconstrained update of the con¬ 
figuration x n . Note that this allows to remove en¬ 
ergy variations due to the discretization of both the 
Hamiltonian and fluctuation/dissipation parts. In 
practice, one has to be careful that internal energies 
should however remain positive. 

Enforcing the total energy conservation allows us 
to obtain a well-defined steady state, for which aver¬ 
age properties can be safely computed. For the equi¬ 
librium simulations presented in the remainder of 
Section HV1 we use the SEM, SSA, SER and Hybrid 
schemes presented in Section Mil but corrected by 
the projection procedure on the internal energies. 
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D. Timestep bias on average properties 
1. Sequential simulations 

The results presented in this section are obtained 
by computing averages over one long trajectory of 
physical time tf = 1000. Figure [4] shows the average 
temperatures obtained from the estimators ([9j, (1101) 
and 111 1 ) for o = 4, while results for the smaller fluc¬ 
tuation strength a — 2 are reported in Figure [5] 



At 


Figure 4. Numerical estimations of the equilibrium tem¬ 
perature as a function of the timestep for a — 4. Top: Ki¬ 
netic temperature ([9} and potential temperature 1101 . Bot¬ 
tom: Internal temperature 1111 



Figure 5. Numerical estimations of the equilibrium tem¬ 
perature as a function of the timestep for a — 2. Top: Ki¬ 
netic temperature l[9j and potential temperature 1101 Bot¬ 
tom: Internal temperature 1111 

As expected, all the estimations of the tempera¬ 


ture extrapolate to the same value as At —» 0, for all 
schemes and for all components of the energy. The 
second point is that the systematic biases depend on 
the timestep At as expected, but also on the value 
of o\ larger values of a lead to much larger errors. 
The third point is that, in opposition to the biases of 
SSA which remain very small even for large values 
of At and a , those of SEM and SER are much larger 
e.g, the bias on the estimation of the kinetic tem¬ 
perature for SER and SEM at At = 0.008 and a = 4 
is around 30% of the extrapolated reference value. 
Moreover, in the case of large At and large a (typ¬ 
ically a — 4 and At > 0.006), SER and SEM biases 
are larger than the linear increase in At observed for 
smaller timesteps. In addition, for internal tempera¬ 
tures estimations, the biases between SER and SEM 
deviate at larger timesteps. 


2. Parallel simulations 

Thermodynamics averages in this section are com¬ 
puted over a single trajectory of total time tf = 450. 
The system was decomposed on N proc = 8 processors. 
Figure [6] presents the temperature estimations as a 
function of the timestep for a = 2, compared to ref¬ 
erence values computed by a sequential SSA simula¬ 
tion of a system of N = 4000 particles. 

Simulation results for a = 4 are quite similar to 
sequential results (although with larger biases), ex¬ 
cept for SER. Indeed, the scheme is not stable even 
for small timesteps since the energy reinjection pro¬ 
cedure may lead to negative internal energies. This 
is related to the lower value of C v used for parallel 
simulations. 

CPU timings are reported in the second column 
of Table [I] For the comparison to be fair, all the tim¬ 
ings have been measured on a single core. We also 
compare these timings to reference timings from a 
sequential SSA simulation (third column of Table Q}. 
The last column of Table Q] gives the required CPU 
time of a simulation of the same system as in Fig¬ 
ure [6j using a timestep chosen in order to obtain a 
given accuracy on the estimation of the kinetic tem¬ 
perature. For each scheme, this timestep, denoted 
Afi; m , is taken such that the bias of the kinetic tem¬ 
perature estimation is equal to e lim , where e lim is de¬ 
fined as the largest bias of the Hybrid kinetic tem¬ 
perature estimations of Figured The biases are com¬ 
puted as e At = |(T kin ) A f - (T kin ) A f =0 |, where <T kin ) A , =0 
is obtained by extrapolating the results of Figure [6] 
to At = 0. The timings of the last column of Ta- 
ble[IJ denoted T c ^, are then obtained by multiply¬ 
ing the timings of the second column of Table[I]with 
the number of iterations necessary to reach tf — 100 
when using the timesteps Aq im . The choice of the 
kinetic temperature is arbitrary (as is the choice of 
e lim). However, timings related to other observables 
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Figure 6. Numerical estimations of the equilibrium prop¬ 
erties as a function of the timestep for a = 2 obtained 
with parallel simulations. Top: Kinetic temperature esti¬ 
mates CD- Middle: Average potential energy (U(q)). Bot¬ 
tom: Internal temperature estimates II It . SSA estimates 
have been computed on a system of N = 4000 particles on 
a single processor. 
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show qualitatively similar results. The conclusion is 
that the Hybrid scheme allows to achieve an accu¬ 
racy comparable to SSA for a given CPU cost, while 
being easily parallelizable. On the other hand, the 
SER scheme on its own, or SEM, are not sufficiently 
accurate to provide interesting alternatives. 


scheme 

^cpu 

ratio to SSA 

T’kin 

L cpu 

SSA 

3.72 x 10“ 5 s 

1 

1.82 s 

SEM 

3.35 x 10 -5 s 

0.90 

12.09 s 

SER 

7.85 x 10 -5 s 

2.11 

28.34 s 

Hybrid 

6.40 x 10“ 5 s 

1.72 

2.13 s 


Table I. Column 2: average CPU time per particle per 
timestep and per processor. Column 3: timings of column 
2 divided by the SSA timing. Column 4: Required CPU 
time corresponding to simulations with timesteps such 
that the kinetic temperature estimation has a given bias. 
This bias is chosen to be the largest bias of the Hybrid ki¬ 
netic temperature estimations of Figure[6] 


We see from Figure \6\ that, as in Section IIV D 11 
SEM and SER estimations are very similar. However, 
Table Q] tells us that SER is much slower than SEM, 
and that using SER takes more than twice as long as 
using SEM in order to reach a given accuracy. This 
subperformance of SER is tempered by the fact that 
we correct here all schemes with the projection of 
Section llV Cl For the Hybrid scheme, Figure[6]shows 
that it yields much smaller biases on average prop¬ 
erties, comparable to those of SSA. This increased 
accuracy more than compensates the extra computa¬ 


tional cost compared to a very simple scheme such 
as SEM, as it can be seen in Table Q] 


E. Influence of the parallelization on the Hybrid 
scheme 

There are no other parameters than the timestep 
for SSA, SER and SEM schemes. This is not the case 
for the Hybrid scheme, where the spatial redistribu¬ 
tion of the simulation box between processors affects 
the computation. It is thus necessary, in order to fur¬ 
ther validate the Hybrid scheme, to study the influ¬ 
ence of the ratio N/N pwc , where N proc is the number 
of processors used for the simulation. Indeed, the 
ratio N/Np IOC directly determines the fraction of in¬ 
teractions treated with SER or SSA. 

In order to study this influence. Hybrid simula¬ 
tions are performed on N proc = 8 processors for sys¬ 
tems in the same thermodynamic state as those of 
Section IIV D 21 The timesteps of the simulations 
range from At - 5 x 10 -4 to At = 2.5 x 10 -3 and the 
number of particle from 40x10x10 to 160x10x10. 
Note that the parallelization is performed in the x di¬ 
rection, and we only change the number of particles 
along this axis. Averages obtained with the Hybrid 
scheme are presented in Figure [7J together with the 
reference SSA and SER estimations of Figured The 
N/Nproc ratios according to the system sizes are dis¬ 
played in Table HI1 



Figure 7. Kinetic temperature estimations for Hybrid sim¬ 
ulations with a number N x of particles along the x axis 
equal to 40, 100 and 160. The SSA and SER estimations 
displayed are those reported in Figure [6] 


system size 

N/N proc 

40x lOx 10 

500 

100 x 10 x 10 

1250 

160 x 10 x 10 

2000 


Table II. N/N proc ratio according to the system size. 
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We notice in Figure [7] that, as expected, Hybrid 
estimations go from those of SER towards those of 
SSA as N/N pioc increases. In the particular situation 
considered here, SER and SSA estimations have bi¬ 
ases of opposite signs that compensate in the Hybrid 
scheme, thus leading to biases smaller than those 
of SSA. Averages of the internal temperature exhibit 
similar behaviors as in Figure [7] SSA and Hybrid 
biases on the potential energy are almost null what¬ 
ever the ratio N/N pTOC . 


V. ACCURACY ON DYNAMICAL PROPERTIES: 
THE EXAMPLE OF SHOCK WAVES 

We study in this section the dynamical properties 
of systems out of equilibrium, using the schemes 
presented in Section [TIT] The system we consider is 
composed of N = 450,000 DPDE particles. In or¬ 
der for the numerical values to be easier to interpret 
from a physical viewpoint, we work in this section 
in physical units. The Lennard-Jones parameters 
in (122b are those of Argon: = 1.65710 -21 J and 

r\\ = 3.40510 -10 m, while the masses of the particles 
are set to m - 6.634 x 10 -26 kg. The initial condition 
is obtained by melting a face centered cubic crys¬ 
tal composed of 250 x 15 x 15 = 112,500 unit cells, 
at density p = 1228 kg.m -3 , with periodic bound¬ 
ary conditions. We use 50 x 2 x 2 = 200 cores for 
the simulations, and various timesteps ranging from 
At = 10 -15 s to At = 4xl0 -15 s (the value At = 10 -15 s 
corresponds to At* = 4.510 -4 in reduced units). A 
timestep of At = 10 -15 is similar to the one used to 
integrate Hamiltonian dynamics of shock waves in 
Argon, but is already one order of magnitude larger 
than for full atom simulations of molecular systems. 
Finally, we set C v = 1.381 x 10 -22 J.K -1 . 

An important point in this section is that we no 
longer correct the schemes by the projection of Sec- 
tion llV Cl since we consider a nonequilibrium system 
which is inhomogeneous in space. This makes in¬ 
deed the global energy reinjection procedure dubi¬ 
ous. 

We equilibrate the system at a given tempera¬ 
ture T re f = 1000 K by superimposing to the DPDE 
equations a Langevin thermostat on the momenta. 
With the notation introduced in j6j, we set the fric¬ 
tion to y — 10 -13 kg.s -1 (which corresponds to a = 
5.255x 10 -17 kg.m.s -372 or a = 7.37 in reduced units). 
Once equilibration is performed, we remove the pe¬ 
riodic boundary conditions in the x direction, and 
put two walls of fixed Lennard-Jones particles of in¬ 
finite masses on the sides of the simulation box to 
confine the system. We next set the system in mo¬ 
tion towards the left wall by adding a velocity equal 
to Up = -2000 m.s -1 to all the particles and to the 
right wall. The velocity of the right wall is main¬ 
tained at lip throughout the simulation. This leads 


to the apparition of a shock wave propagating from 
the left to the right of the simulation box, with an 
average null velocity in the shocked state. 

Figure [8]displays simulation results obtained with 
the SEM scheme with a small timestep At = 10 -15 s. 
Note that no stationary profile is obtained as the in¬ 
ternal temperature increases in time. A similar in¬ 
crease is observed for the kinetic temperature. This 
increase is an artifact of the numerical scheme: we 
indeed confirmed by additional simulations (not re¬ 
ported here) that the increase in the temperature is 
more pronounced for larger timesteps. 



x (m) 


Figure 8. Instantaneous internal temperature profiles on 
the slabs of the cross-section ( y,z ) for the SEM scheme in¬ 
tegrated with a timestep At = 10 -15 s. 

On the other hand, for shock simulations per¬ 
formed with the SER or Hybrid scheme, no such 
drift is observed. Figure [9] presents Hybrid internal 
temperature profiles for various times, for a simula¬ 
tion performed with the same timestep At = 10 -15 s. 
Let us also mention that results obtained with SER 
(not reported here) are completely comparable to the 
ones presented in Figure [9] 



x (m) 


Figure 9. Instantaneous internal temperature profiles on 
the slabs of the cross-section (y,z) for the Hybrid scheme 
integrated with a timestep At = 10 -15 s. 
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Finally, the results are very similar for Hybrid sim¬ 
ulations with timesteps up to At = 4 x 10 _15 s, as re¬ 
ported in Figure flOl 



x (m) 


Figure 10. Instantaneous internal temperature profiles on 
the slabs of the cross-section (y, z) for the Hybrid scheme 
integrated with timesteps ranging from At = 10 -15 s to 
Af = 4x 10 -15 s at time t = 40 ps. 

This shows that the schemes we developed for 
equilibrium simulations, in particular the Hybrid 
scheme, can be used without any projection proce¬ 
dure to simulate dynamical evolutions out of equi¬ 
librium, with timesteps which are even larger than 
the ones used for standard deterministic dynamics. 


VI. CONCLUSION 

We presented two new numerical schemes for the 
integration of DPDE, based on a splitting strategy 
between the Hamiltonian part of the dynamics and 
the fluctuation/dissipation. The first scheme, SER, 
is constructed by decomposing the kinetic energy 
variations into elementary contributions related to 
each particle. This allows to adjust the internal ener¬ 
gies in order for the total energy to be constant. The 
decomposition into elementary contributions can be 
seen as a discrete counterpart to the Ito calculus 
used for the derivation of the dynamics in the inter¬ 
nal energies. The other scheme, namely Hybrid, is a 
blending of SSA and SER scheme, where a SER dis¬ 
cretization is applied to particle couples belonging 
to different processors. 

We carefully compared the ability of these new 
schemes to compute equilibrium properties. All 
schemes exhibit some spurious energy drift over 
time, which can however be straightforwardly cor¬ 
rected by an appropriate projection (performed here 
by adjusting the internal energies^). The scheme 
obtained by applying a projection to Hybrid has an 
accuracy quite comparable to the reference accuracy 
provided by SSA. Therefore, Hybrid is a realistic 
option for a massively parallel DPDE simulation - 


in contrast to other straightforwardly parallelizable 
schemes where the fluctuation/dissipation is inte¬ 
grated with a simple Euler discretization. The SER 
was found to be less satisfactory since the individ¬ 
ual components of the energy show quite large bi¬ 
ases. In addition, SER was found to be less stable 
than SEM or SSA when corrected schemes are consid¬ 
ered. However, its blending into the Hybrid scheme 
cures this somewhat degraded accuracy and stabil¬ 
ity. Moreover, SER, in opposition to SSA or Hybrid, 
is threadable, thus compatible with the future archi¬ 
tectures of supercomputers. 

We finally validated the dynamical behavior of 
our schemes on a non-equilibrium simulation of a 
shocked Lennard-Jones material performed on hun¬ 
dreds of processors, a simulation impossible with a 
naive discretization of DPDE. In our opinion, this 
represents the strongest contribution of this work 
since it shows that the schemes we introduced, in 
particular the Hybrid scheme, allow to easily simu¬ 
late dynamical evolutions of large scale systems on 
massively parallel architectures while at the same 
time being very accurate on the prediction of equi¬ 
librium properties. 
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Appendix A: Stochastic Velocity-Verlet 


The Stochastic Velocity-Verlet algorithm, denoted 
SVV, is an adaption to the DPDE setting of the well- 
known Velocity-Verlet algorithm to integrate Hamil¬ 
tonian dynamics^. It consists in using a Strang split¬ 
ting of DPDE, where the positions are updated 
once with a timestep At, and the momenta p,- and 
internal energies e, are updated twice, both before 
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Note that the same random number is used in the 
first and second updates of the momenta. However, 
the fluctuation/dissipation part has to be computed 
twice per step. 


Appendix B: Consistency of numerical schemes 

There are several notions of consistency for 
Stochastic Differential Equations (SDEs), see for in¬ 
stance Ref. [331 In this work, we consider the con¬ 
sistency of numerical schemes in the sense of weak 
errors. As a corollary, this allows to estimate errors 
on average properties when the schemes are ergodic. 

In order to be more precise, we define the evolu¬ 
tion operator associated with the SDE JTJ: 

(e ,c A){x) = IE[A(x f )|x 0 = x], (Bl) 

where x = (q,p,e) and IE denotes the average over 
all realizations of the Brownian motions in CD- In 


words, the evolution operators give the average 
value at time t of an observable A, for a system 
started in a given configuration x 0 at time 0. The op¬ 
erator C is called the infinitesimal generator, and is 
the adjoint of the Fokker-Planck operator associated 
with the dynamics. For DPDE, it writes^ 


N 

£ =LS' v *- v ^' v p 


i=l 


N 




7=1 
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~7ij v ij + -^-Ai 


2 ! 7 


where 


A 7 = (v p ,-v P; ) + f (d Ei + d £j ). 


(B2) 


(B3) 


The numerical schemes we consider in this work 
can be represented by a mapping ® At (x, G) as x n+1 = 
® At (x", G"). The evolution operator P Af associated to 
a numerical scheme is defined as 


(P At A) (x) = IE [A(x’ I+1 )|x"= x], (B4) 

where IE now is an average over all the random vari¬ 
ables involved in the computation of x n+1 . 

The weak error of a numerical scheme is the max¬ 
imum error between the time averages of the exact 
solution and the ones predicted with the numerical 
scheme. A method is of weak order rj if, for any ob¬ 
servable A and for any simulation time T = Nj ter Af, 
there exists K > 0 such that, for At sufficiently small, 


sup 

0 <n<N it 


E[A(x MAt )]-E[A(x")] <KAt\ (B5) 


where averages are taken for a given initial condi¬ 
tion Xq = x°. In order to be of weak error rj, a numer¬ 
ical scheme should satisfy^ 

{Pm A) (x) = (e £At A)(x) + 0(Af' +1 ). (B6) 

In fact, when this condition holds, it can be proved 
that, provided the numerical scheme and the contin¬ 
uous dynamics are ergodic, the invariant measure of 
the numerical scheme is correct up to errors at most 
0(Af'i) (and sometimes less, see for instance Ref. fd3 
for a discussion in the context of Langevin dynam¬ 
ics). Therefore, the exponent rj in 11131 is determined 
by expansions such as IIB61 . 

We prove below that (1B61 holds with )/ = 1 for SER 
(see Section lB~2ll . This result is obtained by a general 
approximation result for splitting schemes, recalled 
in Section IBTI 


1. A general consistency result 

We recall in this section the general strategy for 
proving dB61l with rj - 1 for splitting schemes. We 
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assume that the dynamics can be decomposed into m 
elementary dynamics. This induces a decomposition 
of the generator as 

C = Ci +... +C m . 

Let be the evolution operator associated with 
the numerical integration of the elementary subdy¬ 
namics with generators If the splitting is done 
according to the Trotter formula^ 

P&t - Pi,At ■■■ (B7) 

and if (fB6t holds for every subdynamics with 
i] = 1, then dB6b holds for the complete numeri¬ 
cal scheme P& t (see for instance the discussion in 
Ref. HI). 

The Hamiltonian part of DPDE is integrated with 
a Velocity-Verlet algorithm 29 , and the corresponding 
evolution operator turns out to be correct at second 
order. Therefore, in order to prove the weak consis¬ 
tency of order 1 for the schemes under consideration, 
it suffices to prove that the fluctuation/dissipation is 
integrated with a scheme of weak order 1. 

2. Consistency of SER 


vanish in average: = 0. The random variables 

B"- involve terms containing products of Gaussian 
variables, but their averages can be easily computed 
as E [G? j GZ l ] = d6 ik 6 jl . A simple computation then 

shows that E^B"J = da 2 x(r"j)/m. 

Let <&*?(*, G) = (<hf’ p (p,G),®*f' £ ( £ ,G)) be the 
result of an Euler-Maruyama discretization of the 
fluctuation/dissipation dynamics (1151 (i.e. the 
scheme U3). The previous computations show that 
SER can be seen as a perturbation of the Euler- 

Maruyama scheme: p' ,+1 = 0^' p (p", G"),-, and 
£ n + l =0 EE, £(£)!) G«). 

- \ ~ At + C«.At 3/2 + 0(At 2 ). 

i =i ' m ' 

i** 

Since the Euler-Maruyama scheme is weakly consis¬ 
tent of order 1, and since the two dominant terms in 
the above expansion vanish in average, a simple Tay¬ 
lor expansion then shows that the integration of the 
fluctuation/dissipation part in SER is also weakly 
consistent of order 1. 


A simple computation shows that (121 ll implies 
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